A scheme for fast exploratory simulation of 

azimuthal asymmetries in Drell-Yan 
experiments at intermediate energies. The 
DY_AB Monte Carlo event generator. 



A. Bianconi 

Dipartimento di Chimica e Fisica per I'Ingegneria e per i Materiali, Universitd di 

Brescia, 1-25123 Brescia, Italy, and 
Istituto Nazionale di Fisica Nucleare, Sezione di Pavia, 1-27100 Pavia, Italy 



Abstract 

In this note I report and discuss the physical scheme and the main approximations 
used by the event generator code DY_AB. This Monte Carlo code is aimed at pre- 
liminary simulation, during the stage of apparatus planning, of Drell-Yan events 
characterized by azimuthal asymmetries, in experiments with moderate center of 
mass energy yfs « 100 GeV. 
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1 Introduction 



Here I discuss the general (physical) scheme of the series of event generators 
DY_AB, concentrating specifically on the last version of this code DY_AB5. 

The event generator DY_AB5 is a generator of dilepton Drell-Yan events [1] in 
hadron-hadron, hadron- nucleus and hadron- (partly polarized molecular tar- 
get) collisions. It is aimed at fast preliminary simulation of that subset of 
Drell-Yan experiments, where 
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(i) the center of mass energy is "intermediate" (from a few to some tenths 
GeV); 



(ii) the projectile may be any light hadronic species (charged pion, proton, 
antiproton), possibly polarized; 

(iii) the target is in general a molecular species, with partial normal polariza- 
tion of some of its component nuclei; 

(iv) the final leptons present azimuthal asymmetries and these asymmetries 
are the goal of the measurement. 

Several experimental proposals have been presented or are in preparation in 



The main difficulty of such experiments is the need to select regions of the 
overall phase space where the event rates are small (in particular: transverse 
momentum over 2 GeV/c), and where two event numbers (e.g.: before/after 
reversing spin) must be compared to identify small asymmetries. It is essential 
to understand from the very beginning which overall event numbers are needed 
to reach a satisfactory population of the interesting subregions. 

The here discussed code is aimed at such "preliminary" investigations, for ex- 
perimental planning only. Since it is based on strong phenomenological com- 
ponents, it is not suitable as it is for theoretical analysis at quark-parton level. 

Up to date, five (private) versions of this code, named DY_AB1, 2, 3, 4, and 
5, have been used and tested by the author and by other users both for 
phenomenological publications [71IS1PllUfl4fllj and for exploratory simulations 
aimed at experimental proposals [21315] (see e.g. [T2], [T3]). 

The latest release DY_AB5 is publicQ. 

This code is not a multi-purpose code. Its main advantages are in its specificity, 
and are: (i) easy insertion of new parametrizations for distribution functions 
associated with azimuthal asymmetries, (ii) easy control and modification of 
the code, (iii) possibility of simultaneous treatment of events in the Collins- 
Soper reference frame and in the fixed target or collider frame, (iv) fast gener- 
ation of events, (v) satisfactory phenomenological reproduction of transverse 
momentum distributions. 

The present note is not the "readme" handbook of the code. The code itself is 
normally accompanied by a readme file supplying user help. Here I discuss the 
physical scheme used for the event generation. This is inspired by the parton 
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model, but some simplifications or phenomenological parameterizations have 
been introduced into the standard relations for the cross section. 

The reasons behind these simplifications are two: 

(i) This code is not aimed at improving the theoretical understanding of quark- 
quark interactions; it is used for reproducing as realistic as possible event 
distributions and associated errors, in measures where some gross features of 
the data are already well known, while other ones are largely unknown. 

(ii) At the present stage, the real point is to understand whether or not certain 
measurements will be possible, or at which extent they will be possible. This 
requires a huge amount of exploratory simulations, to be run in the smallest 
possible time and with the maximum possible flexibility. 

I hope this presentation clarifies in which frameworks the code can be used, 
or should not be used. 



1.1 Development notes 

This code has been written in C++ since the first version. It began as a 
toy model Drell-Yan generator, aimed at fast exploratory simulation for the 
Drell-Yan measurement within the PANDA experiment [12j. After the very 
first applications, the number of options has increased exponentially. 

It was initially used by people of some experimental collaborations, in a form 
that permitted them to handle input in a simple way, assuming that they 
would not need touching the code. This has shown to be unrealistic. On the 
other side, unrealistic as well has been the hope that people could be made 
able to modify the code themselves, without interacting at all with the author. 
And also attempts to organize a "once and for all" form of the code have failed, 
just for the fact that the field is quickly evolving. For example, it is difficult 
to find a "universal" form for the new distribution functions that one could 
like to insert in the next five years. 

So, the general idea is that apart for a central core of classes/functions there is 
nothing sacred in the code structure, and that users should be able to modify 
the code via an as small as possible interaction with the author. 

The first versions like DY_AB1 fully exploited the possibility of writing com- 
plex hierarchies, offered by C++. In DY_AB4 this structured form was aban- 
doned, but for efficiency purposes massive use was made of pointers. DY_AB4 
is well tested, and is the most efficient code of this series. A short presentation 
of this code may be found in [T2] . Its main disadvantage was hard readability, 
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since to increase efficiency it exploited systematically the fortran-style tech- 
nique of organizing big data structures, with functions working on these data 
without explicitly getting them as arguments (the "common" areas in fortran; 
in C++ the same is obtained via pointers to data classes). The absence of 
explicit arguments in the function calls makes the code hierarchy difficult to 
see at first reading. 

DY_AB5 is less efficient (about 30 % more time-consuming). The advantage is 
that it is much easier to read and modify (pointers have mostly disappeared). 
It offers more pre-cooked options as far as unpolarized and single-polarization 
Drell-Yan are concerned. 



2 General theoretical notes 

2.1 General scheme 

This code typical cycle consists of the following steps: 

1) Random event generation in the Collins-Soper frame for a specified class of 
projectile and of target hadrons (e.g. negative pion and polarized proton). 

2) Transformation of these events first to the hadronic center of mass frame 
("collider frame"), and next to the fixed target frame. 

3) Loop over events taking into account the composition of the (molecular) 
target in terms of different nuclei. 

4) If required, a number of repetitions of a multi-event simulation is performed, 
possibly with spin reversal. 

5) If required, statistical analysis of azimuthal asymmetries is performed, cal- 
culating averages and fluctuations of the results obtained in the simulations 
of step (4). 

Here I discuss the general problems and some implementation detail connected 
with steps (1) and (2). 

The exact and complete cross sections on which the event generation is based 
may be found in [H]. There, two simulation schemes are compared. DY_AB5 
is based on what is named "scheme II" in ref. [13] ■ An alternative code has 
been prepared by this author, based on "scheme I" , and in [TJ] the differences 
in the outcome are shown. Also, the reasons are discussed why "scheme I" is 
more proper and supposed to take over in the long run, but not yet well suited 
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given the present state of the art on the phenomenology side. 



Formulas reported in [JJ] and implemented in DY_AB5 are rather long. Here, 
a simplified version of these relations is discussed, to clarify the general cross 
section form, the main approximations contained in it, and the way the code 
exploits it. 

In the language of [H] "Scheme I" is the parton model cross section for Drell- 
Yan events with most of the necessary details: several products of distribution 
functions associated with unpolarized and polarized partons and hadrons, and 
with all quark and antiquark flavors, are summed . 

"Scheme II" exploits some more approximations: (i) A cumulative event dis- 
tribution is factorized out of the sum of all terms. This distribution is built 
by a set of simplified parton distribution functions /, and is phenomenological 
in the transverse momentum dependence, (ii) Generalizing an approximation 
method that may be found e.g. in|15j all the least known distribution functions 
h (those associated with azimuthal asymmetries) are expressed as ratios h/ f. 
(iii) Each such h/ f term is added to the sum, "valence- weighted" by ratios 
like 4/9 etc. 



2.2 The overall cross section 



The relations discussed in the following may be found, e.g., in [17] and [18] . 
The problem is that in these two works several of such relations assume sys- 
tematically different forms. Since this difference is present in all the relevant 
literature, we name these two works "ref.A" and "ref.B" and systematically 
report the differences. 

I name "mass" M the invariant mass of the lepton pair (it is also named Q 
in some references). The indexes "1" and "2" represent the target and beam 
hadrons. 

The Drell-Yan differential cross section can be written in an approximate 
factorized way, inspired by the parton model (see [16], chapter 5, and refs.A 
and B): 



^-■S(r,X F ).S'(P t ).A(e, 



drdXpdPtdn 
or equivalently in the form 
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da 



K(t) 



■ s(x 1 ,x 2 )-s'(p t )-A(e,<p,<i>s). 



dX 1 dX 2 dP t dn 



s 



As customary, s is the squared invariant mass of the two colliding hadrons. 



The scaling assumption means that the only dependence of eq(l) or (2) on s 
should be contained in the 1/s term. 

Xp and t are invariant adimensional variables associated with the beam-axis 
momentum projection, and with the virtuality, of the virtual photon produced 
inp + p— > X + 7* — > X + The pair r, X F can be substituted by 

the equivalent pair X-y, X 2 . Then S and S are related by a Jacobian factor. 
The definition of r, Xp, X\ and X 2 is given below, where it requires some 
discussion. 

Pt is the transverse 2— dimensional component of the 4-momentum of the 
virtual photon, with respect to the beam axis, (it is also named q t , or Q t ). 

The angular variables 6, 0, <p s appearing in A{6, <p, <f> 8 ) are measured in a 
reference frame where the virtual photon is at rest. In this frame 9 and are 
the polar and azimuthal angles of the momentum of one of the two muons, 
while <f) s is the azimuthal angle of the target spin. These variables are discussed 
later, in the paragraph on the angular distributions. 

In the right hand sides of equations (1) and (2) S(...) and S(...) differ because 
of the direct substitution r = t(X±, X 2 ), Xp = Xp(Xi,X 2 ), and because S(...) 
= JS(...), where J is the Jacobian of the coordinate transformation between 
drdXp and dX\dX 2 . 

The coefficient K = K{f) is normally named "K-factor" and is predicted to 
be 1 in the parton model. Actually it is neither 1 nor constant (it is ~ 2). 
Traditionally K contains all the parton model violations, which are so kept 
apart from the rest of the cross section. Summarizing the PQCD corrections 
into a single r— depending factor is at a certain extent justified for r far from 
its kinematic boundaries and 1 (see [TB], subsections 5.2, 5.3, 5.5, and [19j), 
and for moderate transverse momenta << M. At these conditions the parton- 
parton — * 7* + X cross section is dominated by those terms where X subtracts 
no invariant mass from the parton-parton — > 7* transition predicted in the 
plain parton model. 

The above cross section factorization is formally exact within the parton 
model, but in the above reported form it hides that S'(P t ) is (weakly) de- 
pendent on t and Xp, and that A(9, 0, (p s ) depends on Xi, X 2 , P t , M. 



s = 



{E CM f = 2m p (m p + E PtLAB ). 



(3) 
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In the code these dependencies have been taken into account. The previous 
equations are however written in such a way to focus on the assumed variable 
hierarchy that allows for a kinematic separation of the S, S' and A contri- 
butions (see e.g. refs.A and B for examples of the experimental procedure to 
extract these terms from incomplete phase space): 

(i) S(Xp,r) does not depend on P t and on the angular variables, so for any 
assigned Xp,r (or equivalently, Xi,X 2 ) it can be calculated from the cross 
section integrated over all the P t ,9,<p phase space and summed over spin; 

(ii) S'(P t ) = S'(Xf, r, P t ) does not depend on the angular variables, so it can 
be determined by 9, (f> integration plus a sum over spin. 

The function S'(P t ) and A(9, 0, <f> 3 ) are defined so that 

Js\P t ) = 1, jA{9AAs)d^ = 1, (4) 

integrating over all the phase space, for any given <p s . So the total a is just 
the integral of K{t)S{X u X 2 ). 



2.3 The longitudinal term S(Xi,X 2 ): definitions and dangerous ambiguities. 

We can describe the meaning of r and X F by the following relations: 

M 2 




Ml 



(5) 



max 



Irnax ) CM 

yfr is the ratio of the virtual photon virtuality M 2 to its kinematic maximum 
s, reached in an exclusive pp — > l + l~ annihilation into dilepton. Xp is approxi- 
mately the ratio of the beam axis component of the virtual photon momentum 
(in the hadron collision CM) to its kinematic maximum. The precise defini- 
tion of Xp is not univoque in the literature, as discussed below. Whatever the 
exact definition, Xp and r are normally defined as measurable scalar func- 
tions of the projectile and target four-momenta. Alternatively, they can be 
substituted by their combinations Xi, X 2 , whose approximate meaning is the 
ratio of the longitudinal component of each of colliding quarks to the parent 
hadron momentum in the reaction CM: 

X t » f ffi^O • (7) 

V rz)hadron J CM 
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For X\ and X2 several definitions can be found, all approximately equivalent 
at large s and M. These definitions fall into two groups: (a) "theoretical" 
definitions, given in terms of the (unmeasured) light cone momenta of the 
colliding quarks; (b) "experimental" definitions (as in ref.A or B), which ex- 
press X\ and X2 as combinations of the (measured) variables r and Xp. In 
the "theoretical" case, X, is the ratio of the large light-cone component of the 
i— quark momentum to the corresponding component of the momentum of its 
parent hadron. For a rigorous theoretical definition of X\ and X 2 see e.g. |20j . 

In the high energy limit experimental definitions are supposed to reproduce 
the corresponding theoretical ones, so to access approximately to the quark 
momenta. It must be remarked that this is not the situation, in the some 
portions of the kinematic range considered here. 

The definition of r is the same in refs.A and B, and seemingly "r" means the 
same thing in all the literature on the subject: 

r = M 2 /s (refs. A and B) (8) 



On the contrary, for Xp, Xi, and X2 we have non univoque definitions. Ref.A 
uses 



2 P 

X F = (ref. A) (9) 

X F = X 1 -X 2 , r = X,X 2 (ref.A). (10) 



Ref.B uses 



o p 

Xf = * ; (ref. B) (If 

vM 1 -r) 

X F = (A 1 -X 2 )/(l-r), r = X X X 2 (ref. B). (12] 
The definitions of ref.A are easier to use and more common in the literature 



so the code sticks to them. With them it is necessary to take care with the 
kinematic limits: |Xi?| maa: < 1. 

When comparing differential cross sections referred to the variables X\, X 2 or 
Xp, t, jacobian conversion factors are necessary. 

The differential cross sections of equations (1) and (2) enjoy complete scaling 
properties for the pp case, while in the Tip case a mass-dependent [2TJ term is 
introduced by ref.A, it is important for large X^ and considered in the code. 

Full scaling means that (P t , 9, <fi)— integrated cross sections only depend on 
the X— set variables, apart for an s dependence confined to the 1/s term. 
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Compared data of ref.A (250 GeV/c beam energy) and ref.B (125 GeV/c) on 
7i — p DY scattering show approximate scaling. 

K is assumed to be a function of r in ref.A (and in the code), while most 
experiments (see the DY database from[22]) use it as a constant normally ~ 
2. For large r values, the data by ref.A show that this dependence is relevant, 
and the results of the calculation by [19] support this point. The largest part 
of the events concentrate at the lowest part of the involved r— range (wherever 
it begins), and this may make it difficult to scan a large r— range, to establish 
precise dependencies for K on r. 

As remarked in ref.B, the choice of the K value depends on the choice of 
the normalization for the quark distribution functions, which is not univoque. 
Ref.B reports a detailed and systematic discussion of the different normaliza- 
tion methods and of the consequent changes in the values of the distribution 
function parameters and of K. This can be exploited, with the warning of 
using distribution functions according to the notations of ref.B (see below). 

In the default version of the code DY_AB5, S(t,Xf) has been reconstructed 
using the parameterized form given in appendix A of ref.A, together with the 
kinematic definitions and structure functions contained in the main text of 
ref.A. 

This allowed me to fit 7r _ — Tungsten DY differential cross sections reported 
by that experiment at 252 GeV/c. 

When notations of the two works are conformed each other, the distribution 
functions fitted from ref.A allow to reproduce reasonably 7r~— data reported by 
ref.B at 125 GeV/c. To reproduce the p — A DY data of the same experiment, 
proton quark distribution functions must be used as antiproton antiquark 
distribution functions, and then the reproduction is reasonable. 

In both references and everywhere in the literature S(Xi,X 2 ) has the form 

S(X 1 ,X 2 ) = G(X l ,X 2 )E i F(X l )F(X 2 ) (13) 

where G(Xi,X 2 ) is a kinematic factor proportional to (1/A 1 A 2 ) 2 (ref.A) and 
l/XiX 2 (ref.B). In general its exact form depends on notations and changes 
from paper to paper. The exponent "1" or "2" in the 1/X\X 2 factor is very 
important because it indicates whether the distribution functions F(X) must 
be read as F(X) or as XF(X) (see below). 

F and F are linear combinations of the q/q main distribution functions u(X), 
d(X), s(X). For these we have 
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X-u(X) A =u(X) b , X-d(X) A =d(X) B , etc. (14) 



So, e.g., the normalization / XdX(u + d) = 0.34 in ref.A becomes / dX(u + d) 
= 0.34 in ref.B. The a parameter appearing in the typical parameterization 
u(X) = X a {l — X) 13 changes by one unit passing from ref.A to ref.B, and so 
on. 

The important remark is that this ambiguity q(X) «-> Xq(X) is present 
throughout all the literature on the subject, not only in these works. This 
is a very delicate point and must be taken into account whenever new terms 
are added to the code. 



2.4 The P t dependent S'{P t ). 

The traditional parton model literature is built on the collinear approximation. 
So, for the P t — dependent parts one must rely on phenomeno logical fits. 

Experiments A and B did not impose a low-Pi cutoff, with the consequence 
that their small P t data show a completely different qualitative behavior. Mea- 
sured values of the function S'(P t ) can be seen e.g. in ref.A figs. 23 and 25 (n+p 
case), or in ref.B fig. 9 (ir+p and p+p cases). Since azimuthal asymmetries are 
very small for P t < 1 GeV/c, this difference is not relevant for the purposes of 
planning experiments on azimuthal asymmetries in Drell-Yan. For pions the 
default option is the distribution used in ref.A. 

The code however offers a series of alternatives. The class that handles all 
Pt- distributions is PT2. The use of a specific distribution requires subclassing. 

Some options are already present in the code DY_AB5: 

PT2_01d : public PT2 is further subclassed into three possibilities: 

(1) The distribution by Conway et al[17] to reproduce 7r~— tungsten at 250 
GeV/c. 

(2) The one by J.Webb [23] (E866 collaboration) for proton-nucleus at 800 
GeV/c. |hep-ex/030103T1 

(3) The one by Chang [24J (E866 collaboration) for the same measurement of 
J.Webb at the J /if) mass. 

(4) The class PT2_Simple_Asym : public PT2 of the form NP t n /(P? + P 2 ) m . 
This shape is useful for the azimuthal asymmetry terms (see below). 
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The most relevant features of the measured S'(Pt) are (i) a not too strong 
dependence on Xi,X 2 ,s, (ii) for P t < 2 GeV/c the distribution is not steeply 
decreasing (in case ref.A it is increasing up to 0.5 GeV/c); For P t > 2 GeV/c 
it decreases steeply (but with power law, not exponential); (iii) the average 
Pt is near 1 GeV/c, and as well known (see e.g. [16]) it is larger than in 
lepton-induced DIS and in hadron-hadron semi-inclusive meson production. 

In the preliminary simulation of a Drell-Yan experiment on azimuthal asymme- 
tries, a good phenomenological shape of S\P t ) is a key success factor, because 
measured and/or predicted leading-twist azimuthal asymmetries do increase 
at increasing P t , obliging the experiment to select events at a as large as possi- 
ble P t . However, due to the very fast decrease of S'(P t ) for P t > 2 GeV/c, the 
choice of a too large P t cut-off can make a measurement prohibitive because 
of fast falloff of the event rates at large Pt- 

Ref.A reports an explicit parameterization for S'(Pt) (relative to it— induced 
Drell-Yan). Ref.B reports data and some models for the Pt distributions in 
both pp and in irp DY. In the region P t > 3 GeV/c error bars are too big to 
draw any conclusion, but for P t up to 3 GeV/c np and pp P t — distributions 
are very similar. 

For pion and antiproton projectiles I did not modify the parameterization 
inherited from appendix A (pion) of ref.A and from ref.B (p). The pion one is 
not scale- independent. It depends explicitly on M = st, and produces a slow 
increase of the average < Pf > at increasing s and constant r. 

For proton projectiles, the parameterization by J.Webb (23] is probably more 
recommendable. 

2.5 Isospin/ flavor composition. 

Up to a few years ago the models on the functions associated with azimuthal 
asymmetries didn't go to such details like the Z/A composition of the target. 
For this reason, the first code DY_AB1 did not care isosp in/flavor matters. 
After Hermes and Compass results, some flavor and isospin-dependent param- 
eterizations for single-spin asymmetries have appeared (see e.g. !) 25 2(i 27 ). 
obliging the codes since DY_AB2 onwards to take these problems into account. 

On the experimental side the Z/A composition is important for another reason: 
it determines the effective dilution factor of the target polarization. 

DY_AB5 takes into account events coming from separate pieces of a molecular 
target, taking the individual dilution factor of each nucleon into account. So 
to reproduce a NH 3 target with 85 % polarized H, one may arrange code 
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parameters so to require about 4 events on an unpolarized proton or neutron 
and one event on a polarized proton. After this, any specific sorted event, e.g. 
7r_— neutron, will actually need to be translated into dd or a ss event. 

In DY_AB5 this is taken into account by X— averaged weighting factors. The 
criterion and the weights for the relevant cases are discussed in detail in ref.|14j. 
There, also the errors associated with this technique are shown, in "approxi- 
mate vs exact" scatter plots. 

The point is that a sum of the kind 



where the leading term is common, the weights are constant (and of course 
depend on the projectile-target hadron pair for the sorted event), the terms [...] 
are flavor-specific functions where it is not possible to separate the numerator 
from the denominator. 

Weight factors are needed to compensate the fact that e.g. the asymmetries 
associated with dd collisions have scarce relevance in a process like pp or 7i_p, 
while they have more in n + n. As observed in [T4], the discrepancies between 
the "exact" scheme and the constant-weight scheme are relevant when both 
the colliding hadrons have small X. 

In practice, planned experiments will try to stay as far as possible from this 
region, since common belief is that transverse spin effects are small there. In 
addition, although it would be better to sort events according to the former 
scheme, phenomenological parameterizations do normally extract the ratios 
[S asymmetry I Shading] according to the latter scheme. Paradoxically, as remarked 
in |14j . this makes the latter scheme more proper than the former for a simu- 
lation. 

So, although the author of this note has prepared since long a "DY_AB6" code 
working according with a more exact scheme (the one used for the simulations 
of p3]), such a code is unlikely to be useful for still some time. 

2.6 Angular distributions 

Angles 8, and cj) s in the function A(8, <p, cj) s ) are measured in a reference 
frame with origin in the center of mass of the muon couple. In other words, 



flavor 



E 




(15) 



where each term is flavor-specific, is approximated in the form 




(16) 
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the origin of this frame coincides with the virtual photon position. The axes of 
this frame can be oriented in several ways. One leads to the so-called Collins- 
Soper frame [28J (CS), with z axis parallel to the difference of the momenta 
of the projectile and of the target nucleon. The transverse axes are oriented 
so that the xz plane contains the virtual photon momentum. Other common 
alternatives put the z axis along the beam or target direction in the lepton 
CM. 

The DY_AB5 code sorts events in the CS frame. This fact becomes relevant 
when events are transformed to the fixed target or to the collider frame (see 
related section below). 

In the CS frame, 9 and are the angles formed by one of the muons from 
now onwards), and 4> s is the target spin orientation. 

For a qualitative understanding of the kinematics for most (not all) events, 
one may imagine a virtual photon with transverse momentum not much larger 
than 1 GeV, and \Xp\ not too close to zero, so that the longitudinal momentum 
of the photon is larger than the transverse momentum. Then the CS z axis is 
roughly parallel to the collider one. Also the CS xy plane is roughly parallel 
to the collider xy plane, but the CS x and y axes are randomly rotated by an 
angle (p co u with respect to their collider configuration. This (p co u is the angle 
between the xy— component of the virtual photon and the x axis in the lab. 
The x axis of the CS frame coincides physically with the transverse component 
of the virtual photon 3-momentum. Then the transverse proton spin, which is 
fixed along the x axis in the collider frame, lies on the CS xy plane with angle 

4>s = -<Pcoll- 

These approximations are not used in the code, but can be useful to un- 
derstand the meaning of the employed angles, and to have an idea of the 
distribution of useful events in a collider frame. 

In the Collins-Soper frame the angles 9 and are polar and azimuthal angles 
of the momentum of one of the two leptons and are randomly distributed. 
Also the spin angle <fr s is randomly distributed in the CS frame. In absence 
of information on two of these angles, the third one is homogeneously dis- 
tributed over all of its phase space. However, altogether their distributions are 
correlated by A(9, <p, <p s ) in the cross section. 

In the simulation events, are initially flat-sorted with respect to all the kine- 
matic variables X 1; X 2 , Pt, 9, <fi, (fi s . The sorted events are accepted/rejected 
according to the cross section expressed by eq.(l) where, at the level of single 
spin experiments, 
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A{e,(j),(f> s ) = 1 + cos 2 (6) + K x ' 2 2 ' t} sin 2 {6)cos(2(f>) + 

+ |£ 2 |fl(Xi,X 2 ,P t ,M,0,0,&) (17) 

The unpolarized asymmetry measured in ref.A is contained in the v... term. 
Single spin asymmetries arise from the | S'2 1 -B ( . . . ) term. Two origins for such 
terms have been considered here, corresponding to different origins for the 
azimuthal term. 

The Siversp9] asymmetry considers a term of the form 

B = 2^sin{2e)sin{4> - <p s )H a (X u X 2 ) (18) 
while the Boer-Mulders [T51I30] asymmetry is of the form 

B = -\ s p^sin 2 (e)sin(<P + <p s )H b (X 1 ,X 2 ) (19) 

Any of the above azimuthal asymmetries (unpolarized, Sivers, BMT) can be 
disentangled from the other two by a suitably weighted integration. The 
"statistical analysis" option of the code offers this possibility. 

2. 1 Parameterizations for the nonstandard distribution functions 

As above written, the leading distribution functions are bypassed by assuming 
a phenomenological (correctly behaving and scaling in a wide range of kine- 
matics) form for the "standard" event distribution, i.e. for the (9, 0)— averaged 
part of the event distribution. 

For the "nonstandard" terms (i.e. those ones that produce nontrivial angular 
distributions in the CS frame) the code includes as a first possibility some phe- 
nomenological parameterizations that have become recently available in the 
literature, and as a second option the possibility to chose freely the parameters 
for simple pre-determined shapes. 

The functions associated with the nonstandard terms have been put in places 
where it is easy to find them (in the file c_dy_master.cpp). So, if one wants 
to change radically the shape of these functions it is possible to do it. The 
code has been written assuming that any potential user could be interested in 
adding supplementary terms to this set. In other words, I assume that in the 
next ten years there will be no special reasons to change the "standard" part 
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of the cross section, while updates and news on the asymmetry side will be 
frequent. 

Although it is possible to add new parameterizations, respecting some formats 
makes things easier. Here I would like to discuss these formats. 

1) Let us write the A term in eq.(17) in the most general form 

A( ) = 1 + cos\6) + ]T X 2 ,P t )F'(0)F'U0 s )) (20) 

It must be reminded that according with the scheme adopted in this code 
each F term is the ratio between a term associated with a given kind of 
angular asymmetry and the "standard" term, that carries with itself the 1 + 
cos 2 (9) angular dependence. Since most of the available parameterizations give 
directly such a ratio, this is the most convenient form, as earlier discussed in 
this work. 

2) The code assumes factorization between longitudinal and transverse de- 
grees of freedom, and between terms coming from projectile and target, in the 
functions terms v and B in eq.(17). These functions have the form 



F(X u X 2 ,P t ) = f 1 (X 1 )f 2 (X 2 )f t (P t ), (21) 
3) For the longitudinal components fi(Xi) the code assumes the form 



fi(X) ee NX a (l-Xf (22) 

To insert completely different functions is possible but it is of course much 
easier to change three parameters for any quark flavor. 

4) The transverse term is 



j,^ p ^ = J ' d 2 k 1 d 2 k 2 [g 1 (k 1 )g 2 (k 2 )} asymmetry S 2 (P t - h - k 2 ). ^ 
/ d 2 k 1 d 2 k 2 [g 1 (k 1 )g 2 (k 2 )] standard S 2 (P t - k x - k 2 ). 

where as previously reminded, the numerator is the really asymmetric term, 
the denominator is the leading standard P t — dependence. 

5) For the transverse part it is responsibility of the user to introduce directly 
the convolution of the two transverse momentum distributions corresponding 
to the colliding partons. In other words, the code assumes that one directly 
introduces ft(Pt) hi eq.(20). 
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The last constraint is not a true constraint, since all the parameter izations 
known to me employ gaussian shapes, whose convolution may be easily com- 
puted in analytical way. On the other side, this means spared execution time 
for the code. 

7) Frequently, the result of eq.(23) ft{Pt) "will assume the form that I name 

fsimplei^Pt) • 

fsim P le(Pt) = NP t n /(P t 2 + P*r. (24) 

Because of the common use of Gaussian distributions, this form for the ratio 
f{Pt) is quite common. So the code gives, among other options, this one. 

To implement the calculated Pt— convolutions, or to insert phenomenological 
ones (including the "standard" term S(P t ) in eq.l) the code offers a class PT2 
that may be subclassed two ways: 

(i) Exploiting the class PT_Simple_Asym : public PT2 to directly insert the 
parameters N, P a , n, m, into the form NP™ /{Pf + P 2 ) 171 . 

(ii) Subclassing the class PT2 by another user-assumed distribution. As above 
discussed, the code itself offers several examples of this procedure, i.e. all the 
relevant alternatives offered for the P^-dependence of the leading "standard" 
term. 

The already present default parameterizations correspond, for the Boer-Mulders 
effect, to the x— independent v{k t ) function given in [T5], for the Sivers effect 
to the two alternatives PI2"5] . 

For the k t — unintegrated Transversity distribution the default (N, pa- 
rameters are (1,0,0) and no predefined set was considered. 



3 Interesting Plots 

In the following, some collections of simulated Drell-Yan events are compared 
with the measured distributions, for negative pions on Tungsten. At the kine- 
matics of interest for DY_AB5 the two experiments with far the best statistics 
in the mass range 4-9 GeV/c 2 are E615 at Fermilab and NA10 at CERN. The 
code DY_AB has been based on the scheme and parameters given by E615 
in [IT], however it works equivalently with NA10 data [311132"] . coming from 
nearby kinematic regions. 
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3.1 Comparison with experimental data: E615 

To produce fig.l, DY_AB5 has sorted 100,000 Drell-Yan events for negative 
pions with beam energy 252 GeV on a Tungsten target, in the dilepton mass 
range 4-7 GeV/c 2 . From this set, I extract the subset of events with y/r in 
the range 0.254-0.277 (mass from about 5.5 to 6 GeV/c 2 ). The distribution 
of these events with respect to xf may be compared with the cross sections 
given in [T7], table VI. 



Simulated events vs measured cross section 
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Fig. 1. Filled squares: E615 data for yjr in the range 0.254-0.277. Empty squares 
with joining line: Simulation (see text). 

In fig.l the empty squares joint by a line are the event numbers sorted by 
DY_AB5. Horizontal bars represent the size Axp = 0.1 of each bin, vertical 
bars the statistical error y/N. The full black squares with no horizontal error 
bar are the numbers from table VI of [T7|, rescaled by a common constant 
factor so to transform cross section values into expected values for the bin 
populations. 

Sets of data corresponding to smaller sfr cover a slightly smaller range in 
xp, while at increasing yfr the event numbers filling each experimental bin 
decrease, and error bars increase. For comparison, in fig.2 we report data from 
table VI of p2| in a mass range near the upper edge 9 GeV/c: yfr ranges from 
0.392 to 0.415 (mass between about 8.5 and 9 GeV/c 2 ). Clearly, the error bars 
are much bigger than in the case of fig.l. The corresponding simulated event 
numbers are extracted from a set of 100,000 sorted events between mass values 
7 and 9.2 GeV/c 2 . 
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Fig. 2. Filled squares: E615 data (252 GeV pion beam) for yfr in the range 
0.392-0.415. Empty squares with joining line: Simulation (see text). 



Simulated events vs measured cross section 
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Fig. 3. Filled squares: E615 data for xp in the range 0-0.1. Empty squares with 
joining line: Simulation (see text). 

Data from [17] do not cover xp < —0.1 or —0.2. This is a typical situation in 
fixed target experiments, where x pro j ect u e ~ 1 and x ta rget ~ 0. 

The data and simulations in figs.l and 2 are integrated with respect to the 
transverse momentum of the dilepton pair. In figs. 3 and 4 I report distribu- 
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Simulated events vs measured cross section 



1000- 



200 



800- 



600- 



400- 



n I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i ryj-a^-f tH* 

0.5 1 1.5 2 2.5 3 3.5 





Fig. 4. Filled squares: E615 data for xp in the range 0.6-0.7. Empty squares with 
joining line: Simulation (see text). 

tions with respect to P t , for two assigned ranges of Xp\ 0-0.1 (fig. 3) and 0.6-0.7 
(fig.4). Here data and simulated distributions are integrated with respect to 
the mass (equivalently, to r) over the mass range 4-9 GeV/c 2 . The meaning of 
open and filled squares is the same as in figs.l and 2. Data points come from 
the same experiment of [T7] , (they are reported explicitly in [22] ; hi [T7j figures 
on the P t — distributions are present, but a table of values is not reported) 

3.2 Comparison with experimental data: NA10 

The two richest collections of events at the kinematics interesting here have 
been provided by the collaboration NA10, with negative pions of 194 and 286 
geV on Tungsten [3"T1I3"2"] . For the beam at 194 GeV, the data may be found in 
the final table of [31], while the data relative to the upper energy beam have 
been taken from [22], to which they have been sent as a private communica- 
tion. This experiment did not publish transverse momentum distributions. In 
addition, the covered xf range tends to become rather narrow near the lower 
dilepton mass value 4 GeV/c 2 , where most events concentrate. So the most 
interesting distributions are at larger mass values with respect to E615. 

Data in fig. 5 come from NA10 (194 GeV beam), and the simulated events are 
a subset of 100,000 events sorted by DY_AB5 in the mass range 4-9 GeV/c 2 , 
assuming a negative pion beam of 194 GeV hitting Tungsten. The data in figs. 
6, 7 and 8 all refer to the upper energy beam of NA10, and their simulated 
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Simulated events vs measured cross section 
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Fig. 5. Filled squares: NA10 data at 194 GeV, for y/r in the range 0.33-0.36. Empty 
squares with joining line: Simulation (see text). 
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Fig. 6. Filled squares: NA10 data at 286 GeV, for y/r in the range 0.33-0.36. Empty 
squares with joining line: Simulation (see text). 

counterparts have been extracted from a set of 100,000 events sorted between 
4 and 7 GeV/c 2 , of 100,000 events between 7 and 9 GeV/c 2 , and of 50,000 
events between 11 and 15 GeV/c 2 . 

For the lower energy beam (pions of 194 GeV, meaning s 364 GeV 2 ) fig. 5 
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Simulated events vs measured cross section 
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Fig. 7. Filled squares: NA10 data at 286 GeV, for in the ranges 0.21-0.24 (top), 
0.24-0.27 (middle), 0.27-0.3 (bottom). Empty squares with joining line: Simulation 
(see text). 

Simulated events vs measured cross section 
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Fig. 8. Filled squares: NA10 data at 286 GeV, for sfr in the ranges 0.51-0.54 (top) 
and 0.54-0.63 (bottom). These are equivalent to dilepton masses over the bottonium 
region. Empty squares with joining line: Simulation (see text). 

reports the xp distribution for y/r in the range 0.33-0.36, equivalent to dilepton 
mass between about 6.3 and 6.9 GeV/c 2 . 
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Fig. 6 considers the same \fr range, for the case of the upper energy beam (286 
GeV, meaning s ~ 537 GeV 2 ). In this case this corresponds to dilepton mass 
between 7.6 and 8.3 GeV/c 2 . According with the scaling hyppothesis, the two 
data distributions should be the very similar in the xp range 0-0.6 covered by 
both measurements. 

For \fr < 0.3, the x F range covered by NA10 becomes small and data distri- 
butions are rather flat. Fig. 7 reports data and simulations for the three y/r 
ranges 0.21-0.24, 0.24-0.27, 0.27-0.3, corresponding to masses ranging from 4.9 
to 6.9 GeV/c 2 . 

The compared examination of the three sets of fig. 7 suggests that (i) the 
simulation is worse for smaller \fr, (ii) the K— factor extracted from E615 is 
slightly less steep (in its dependence on r) than the one extracted from NA10. 

The former fact depends on the increasing relevance of sea values for x pro j ec tu e 
at small masses. The default distribution for x w in DY_AB5 comes from E615 
where sea (anti) quarks of the pion play a minor role in the fit of the pion 
distributions. For this reason, data from both NA10 and E615 are reasonably 
fitted for xf > —0.1 with the exception of small-yjr distributions. 

A signal of the relevance of sea partons from the pion side is the shift of 
the xp— distribution peak towards xp = 0. Valence-dominated measurements 
show this peak at xp = 0.1-0.3. When the sea of the pion becomes relevant, it 
is probably better to substitute the default sea pion distribution of DY_AB5 
with more recent ones. Fig. 7 suggests that this may be the case for \fr below 
0.25. 

An interesting feature of NA10 is the presence of a conspicuous set of data at 
masses over the bottonium mass. This is not the situation for which DY_AB5 
has been thought, however it is interesting to try and simulate these events. 
Fig. 8 reports data and simulation for \fr in the ranges 0.51-0.54 (mass between 
11.8 to 12.5 GeV/c 2 ) and 0.54-0.63 (mass between 12.5 and 14.6 GeV/c 2 ). 

We see that DY_AB5 has difficulties in reproducing the shape of these event 
distributions for xf < 0.2. Actually, the "gap" of these data distributions at 
xp ?s looks a little unnatural. More in general, in all the previous figures 
the agreement between montecarlo and data is worse for negative xp, where 
data distributions fall rather steeply. This could be related with the fact that 
negative xp data are at the border of the regions of good acceptance for fixed 
target experiments. 
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3. 3 Sivers asymmetry plots in different calculation schemes 

As observed in Section 2.5, the event cross section may be written in the 
form a = o- Q (X 1 ,X 2 , Pt, 6) [1 + A(X\, X 2 , Pt, 0, ..)], where the former term ex- 
presses that part of the cross section that does not contain azimuthal and spin 
asymmetries, while the asymmetries themselves are contained in A, and where 
A may be approximated in different ways. In particular, in DY_AB5 A con- 
tains flavor weight factors, that are absent a more recent version DY_AB6. In 
DY_AB6 the full cross section o is a sum of independent terms each referring 
to a sea or valence parton. Each flavor contribution carries "its own" asymme- 
try. In DY_AB5 a is a sum of flavor contributions, and A is an independent 
sum of flavor contributions (see[14j), weighted with good sense. 
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Fig. 9. Scatter plot of asymmetries calculated by DY_AB5 (horizontal) and by 
DY_AB6 (vertical), for 7000 sorted n'p Drell-Yan events at s = 100 GeV 2 , in the 
mass range 4-9 GeV/c 2 . From |14j . 

As far as azimuthal/spin asymmetries are neglected, the two codes produce the 
same results (<7o is the same in both cases). When asymmetries are included, 
the scheme implemented by DY_AB6 is more proper. 

Exploiting the equality of do in the two cases, in [H] Drell-Yan events have 
been sorted according with do only, and the corresponding Sivers asymmetry 
(deriving from the A term) has been calculated within the relations of DY_AB5 
and of DY_AB6. The the scatter plot of the two asymmetry calculations for 
each event are reported in several figures. From that work two figures are 
borrowed, showing the two most "extreme" cases. 

Fig. 9 reports 7000 events for n~p Drell-Yan at s = 100 GeV 2 , lepton invariant 
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Fig. 10. Scatter plot of asymmetries calculated by DY_AB5 (horizontal) and by 
DY_AB6 (vertical), for 20000 sorted vr+p Drell-Yan events at s = 100 GeV 2 , in the 
mass range 1.5-2.5 GeV/c 2 . From |14j . 

mass in the range 4-6 GeV/c 2 , transverse momentum in the range 1-3 GeV/c. 
The Sivers asymmetry is parameterized according to [9]. 

Fig. 10 reports 20000 events for n + p Drell-Yan at s = 100 GeV 2 , lepton invari- 
ant mass in the range 1.5-2.5 GeV/c 2 , transverse momentum in the range 1-3 
GeV/c. The Sivers asymmetry is parameterized according to [25] . 

In the former case DY_AB5 and DY_AB6 would produce the same Sivers 
asymmetry. In the latter the difference is striking. As deeply discussed in [T4] 
the difference between the two is proportional to the role of sea (anti) quarks. 
Positive pions on protons, and small dilepton mass, enhance the role of sea 
partons. 

For these reasons, the code DY_AB6 has been prepared, aiming at situations 
where sea partons could become more and more relevant. 

However, the use of DY_AB6 has been up to now restrained by the fact that 
the available parameterizations of functions like the Sivers one are normally 
produced by fitting data within the same "ideological" scheme of DY_AB5. In 
these cases, the use of DY_AB6 would increase errors, instead of decreasing 
them. This is discussed in [14J with plenty of details. 

The main point is that when a "global fit" is undertaken using a scheme like 
the one of DY_AB5, the effect of sea partons is effectively included inside 
valence quark distributions. So, when the results from such a global fit are 
emploied in DY_AB6 for predicting some result, it is very dangerous to add 
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separate sea quark contributions that are already present in valence quark 
distributions. 

The place where DY_AB6 can be more proper is the one of theoretical models 
for the Sivers function, built according with the scheme where each flavor is 
individually considered. But for the aim of modelling an experiment apparatus, 
one normally prefers using phenomenological parametrizations, rather than 
theoretical models. 
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